#!/usr/bin/R

#function which specifies weir equations
###############################################################################################
find.q = function(h.up,h.dn,b) {
	my.q=0
        if (h.up>box[b,box.crest] & h.dn<=box[b,box.crest]) {
                #forward flow, not submerged
                my.q=3.33*box[b,box.width]*(h.up-box[b,box.crest])^(3/2)
                #print(paste(b,'forward',my.q[b]))
        }
        if (h.up<=box[b,box.crest] & h.dn>box[b,box.crest]) {
                #backward flow, not submerged
                my.q=(-1)*3.33*box[b,box.width]*(h.dn-box[b,box.crest])^(3/2)
                #print(paste(b,'backward',my.q[b]))
        }
        if (h.up>box[b,box.crest] & h.dn>box[b,box.crest]){
                #submerged sharp-crested weir formula
                my.q=abs(3.33*box[b,box.width]*sqrt(abs(h.up-h.dn))*
                        ((h.up-box[b,box.crest])+(h.dn-box[b,box.crest])*.381))
                #adjust sign for direction
                if (h.up<h.dn) {my.q=my.q*(-1)}
                #print(paste(b,'submerged',my.q[b]))
        }
        return(my.q)
}

#function that determines flow rates across 'boxs' (weirs/pumps)
###############################################################################################
find.flow = function(my.h) { 
		
        #initialize local variables
	q.old=q.new=0
        my.inflow=array(0,dim=length(my.h))
        my.outflow=array(0,dim=length(my.h))
	tol=1e-2

        for (b in 1:n.boxs) {
		
		upstream=box[b,box.upstream]
		downstream=box[b,box.downstream]
		up.dn = c(upstream,downstream)
                len=box[b,c(box.l.up,box.l.dn)]
                friction=box[b,c(box.f.up,box.f.dn)]
                diam=box[b,c(box.d.up,box.d.dn)]
                gate=box[b,c(box.gate.up,box.gate.dn)] 
		
		if (box[b,box.type]==2) { #box is a pump
                        q.new=box[b,box.pump]
                }
                
		if (box[b,box.type]==1) { #box is a weir
			
			#first, compute flowrate given heads upstream/downstream
			q.new=find.q(my.h[upstream],my.h[downstream],b)

			#now find head loss in pipes, recalc flowrate, iterate until convergence
			i=0
			while(abs(abs(q.new)-abs(q.old))>tol | q.new!=0) {
				i=i+1
				q.old=q.new
				h.friction = 8*friction*len*q.old^2/(diam^5*pi^2*grav)
				Kgate = exp(-7.25*gate+4.55)/60
				h.gate = 8*Kgate*q.old^2/(diam^4*pi^2*grav)
				h.pipe = my.h[up.dn]-h.friction-h.gate #h.pipe=head just before/after weir crest
				h.pipe[h.pipe<0]=0
				q.new=find.q(h.pipe[1],h.pipe[2],b)
				
				print(paste(t,i,h.pipe,q.old,q.new))
				if (i==100) stop()
			}
                }

		my.outflow[upstream]=my.outflow[upstream]+q.new
                my.inflow[downstream]=my.inflow[downstream]+q.new
        }
        return(list(my.inflow,my.outflow))
} #end function find.flow()
###############################################################################################


